//
// Ephi - simulation of magnetic fields and particles
// Copyright (C) 2007, 2008 Indrek Mandre <indrek(at)mare.ee>
// For more information please see http://www.mare.ee/indrek/
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
//

#ifndef __pmap_hpp__
#define __pmap_hpp__

class StaticPMap;
class CubicScalar;

class PMap
{
public:
  PMap (const vect3d& dim, prec_t resolution, prec_t contact = 0.0, prec_t smooth = 0.0, size_t smooth_count = 64);
  ~PMap ();

  bool addQ (const vect3d& pos, const vect3d& v, prec_t q);
  void addUniformSphere (const vect3d& pos, prec_t r, prec_t q);

  void recalc ();

  prec_t getPotential (const vect3d& pos);
  void addEField (const vect3d& pos, vect3d& ef);

  bool inContact (const vect3d& pos);

  void addCubicValues (const vect3d& pos, CubicScalar& epot, CubicVector& cbf, CubicVector& cef);

  prec_t getCrudePotential (const vect3d& pos);
  void addCrudeEField (const vect3d& pos, vect3d& ef);

  prec_t getChargeDensity (const vect3d& pos);

private:
  vect3d bs, hbs, rbs;
  size_t xc, yc, zc;
  vect3d corner;
  prec_t cmag;
  prec_t acc;
  prec_t smooth, smag;
  size_t smooth_count;
  prec_t fsmooth; // fraction of smooth
  size_t xs, ys, zs, xss, yss, zss;
  prec_t sv[8][2];

  // M
  prec_t f35_8R3;
  prec_t f21_4R5;
  prec_t f15_8R7;
  // N
  prec_t f21_2R5;
  prec_t f15_2R7;
  // K
  prec_t f15_R7;
  // pot
  prec_t f35_16R;
  prec_t f35_16R3;
  prec_t f21_16R5;
  prec_t f5_16R7;

  struct particle_t
  {
    vect3d pos;
    vect3d v;
    prec_t q;
    bool dobf;
  };

  struct pnode_t
  {
    particle_t p;
    struct pnode_t *next;
  };

  struct cell_t
  {
    size_t index;
    vect3d mid;
    vect3d qmid;
    prec_t qsum;
    vect3d j;
    vect3d jmid;
    particle_t *particles;
    size_t pcount;
    pnode_t *pnodes;
    struct cell_t *next;
    bool dobf;

    void add (pnode_t *p) { p->next = pnodes; pnodes = p; }
  };

  struct sphere_t
  {
    vect3d pos;
    prec_t r;
    prec_t rr;
    prec_t q;
    prec_t q_by_r;
    prec_t q_by_r3;
    prec_t rho;
  };

  allocator<pnode_t, 1024> a_pnodes;
  allocator<cell_t, 1024> a_cells;
  allocator<particle_t, 1024> a_particles;

  std::vector<cell_t*> hash;
  std::vector<sphere_t> spheres;

  cell_t * findCell (size_t index);
  inline size_t makeIndex (size_t xi, size_t yi, size_t zi) { return zi + yi * zc + xi * yc * zc; }
  void getPos (const vect3d& pos, size_t& xi, size_t& yi, size_t& zi);

  std::vector<cell_t*> cells;

  friend class StaticPMap;
};

class StaticPMap : public StaticElement
{
public:
  StaticPMap (PMap& pmap, prec_t Emul = 1.0);

  void addFields (const vect3d& pos, vect3d& bfield, vect3d& efield);
  void addPotential (const vect3d& pos, prec_t& potential);
  void draw(Screen& screen, const Screen::color& c = Screen::BLACK);
  std::string get_povray_script ();
  bool inContact (const vect3d& pos);

private:
  PMap& pmap;
  prec_t Emul;
};

class StaticPMapCrudeEfield : public StaticElement
{
public:
  StaticPMapCrudeEfield (PMap& pmap);

  void addFields (const vect3d& pos, vect3d& bfield, vect3d& efield);
  void addPotential (const vect3d& pos, prec_t& potential);
  void draw(Screen& screen, const Screen::color& c = Screen::BLACK);
  std::string get_povray_script ();
  bool inContact (const vect3d& pos);

private:
  PMap& pmap;
};

#endif

